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A new approach to the averaged two-particle distribution 
function of a crystalline phase is presented. It includes an 
indirect check of the merit of the Gaussian approximation for 
the local density and a new way to inferring values of the 
thermodynamic variables from simulation data. The equa- 
tion of state and the compressibility of the hard-sphere FCC 
crystal is computed from Tarazona free energy density func- 
tional [Phys. Rev. A 31, 2672 (1985)]. They are in excellent 
agreement with simulation results over the physical range of 
densities up to almost close packing. We also include the 
comparison with the results obtained by two other functional 
approaches which are also excellent. 

PACS numbers: 64.10. +h,64.30.+t,61.66.-f,05.70.Ce 



I. INTRODUCTION 

The density-functional theory applied to non-uniform 
classical fluids has been able to depict a wide range of 
physical properties of simple solid systems. Initially, it 
was intended for developing a theoretical approach that 
was able to describe both fluid and crystalline phases con- 
sistently. Then, it was extended to other crystal prop- 
erties like elasticity and phonon dispersion and to more 
complex systems such as surfaces and liquid crystals (see 
reference [jlj for a recent review on this matter). 

It is well known the crucial role of hard-spheres (HS) as 
the usual reference system of more realistic systems that 
include attractive interactions. Accordingly, a consider- 
able effort has been done to develop free energy function- 
als describing non-uniform HS systems and, at present, 
there exist quite good functional approaches for these 
systems 0. Recently, Denton et al. [|| and Tejero et al. 
j| have analyzed the equation of state of a HS crystal 
obtained from the modified weighted-density approxima- 
tion (MWDA) H and from the generalized effective liq- 
uid approximation (GELA) || respectively. Their anal- 
ysis include densities well inside the stable solid phase. 
This is specially relevant in connection with the solid- 
solid transition recently reported by Bolhius and Frenkel 
|^| in systems of HS with a short ranged attractive inter- 
action (we have learnt that Stell and collaborators had 
already predicted these kind of transitions in the seven- 
ties 0). It happens that, for sufficiently short-ranged 
attractions, this transition occurs at very high densities, 
near close packing. To describe this kind of phenom- 
ena, the theoretical approach for the reference HS system 



should give a reasonable equation of state over all the 
physical density range, even proximal to close packing. 
Thus, in this paper, we evaluate the equation of state of 
a face-centered-cubic (FCC) HS crystal at densities up to 
almost close packing using the Tarazona functional ap- 
proach What it is more important is to compare 
the functional predictions with simulations results. To do 
it, we develop a new method to obtain thermodynamic 
information from simulation data of g(r), the average 
of the two-particle distribution function. This method 
also allows us to infer the ideal and the excess contribu- 
tions (as defined in section II) to the equation of state. 
It also corroborates the merit of the Gaussian descrip- 
tion of the one-particle distribution function. Finally, 
but no less important, the method suggest an interesting 
discussion on the averaged correlation between the par- 
ticles beyond nearest-neighbours. The agreement of the 
Tarazona functional predictions with simulation data is 
excellent up to almost close packing. The same can be 
said for the predictions from MWDA and GELA at least 
up to the densities reported. To test this accordance, we 
compute also the compressibility of the HS crystal and 
compare it with that of simulation, finding an excellent 
agreement. Moreover, we also work out the ideal and ex- 
cess contributions to the pressure and the compressibility 
and compare each one with those corresponding to sim- 
ulation data. Again, the accordance is quite good. Note 
that the evaluation of the equation of state by MWDA 
have some problems at high densities as pointed out by 
Tejero et al. §. 

In the next section, we briefly resume the Tarazona 
functional (TF) together with MWDA and GELA. In 
section III we present the mentioned discussion of g(r). 
The results and a discussion of them are presented in 
section IV. The conclusions are exposed in section V. 



II. FUNCTIONAL APPROXIMATIONS 

Density-functional theories |jj are based on a varia- 
tional principle pT| which allows to propose approxima- 
tions to the Helmholtz free energy. The variational prin- 
ciple establishes that for a given interacting potential, 
fixed external potential and mean density, the Helmholtz 
free energy F[p(r)], as a functional of the one-particle 
density p(r), has a minimum value at the equilibrium 
density. The free energy functional is usually written as: 



F[p(r)} = F id [p(r)]+F ex [p(r)}, 



(1) 



where Fid is the ideal contribution to the Helmholtz free 
energy of a non-uniform system. Its functional form is 
exactly given by: 



pF id [p(r)} = / eMr)(Zn(A 3 p(r)) - 1), 



(2) 



where A is the thermal de Broglie wavelength and (3 = 
1/UbT. The second term of ([!]), F ex , called free energy 
excess, arises from the interacting potential between par- 
ticles. Several accurate approximations have been pro- 
posed for the free energy excess. These are based on a 
mapping of the thermodynamic properties of the non- 
uniform systems onto those of a uniform fluid at some 
effective density (weighted-density). Tarazona ||,|9| pro- 
poses the following expression for the free energy excess 
of the HS solid: 



F ex [p(r)} = / drp(r)A^ ex (p(r)), 



(3) 



where A^ ex (p) is the free energy excess per particle of 
the uniform system at density p and the weighted density, 
p, is given by 



p(r) 



J dr' p{y')w{\t 



r'kp(r)). 



(4) 



The function u>(r), in the integral equation that defines 
p, is specified by requiring that the direct correlation 
function c(r) obtained from the free energy functional 
matches that of the HS liquid in the uniform limit. In 
the MWDA, Denton and Ashcroft use the same mapping 
idea. However, they propose a global map of the free 
energy excess onto the free energy of a unique uniform 
fluid ~ 



(5) 



\p{v)\ = NA* ex (p), 
where the weighted density p is given by 



J drp(v) J dr'p(r')w(\r-r'\;p). 



(6) 



In the GELA, Lustko and Baus use again a global map 
but the weighted density (or effective density) , p, is given 
by the structural mapping p|: 



drp(r) / dr' p(r')c(\r 



drp{r) I dr'p(r')c(r, r'; p(r)), 



(7) 



where c is the direct correlation function. 

An important advantage of the two latter approaches 
is that they require less computational effort. However, 



only the TF is a true functional approach in the sense 
that it is not restricted to macroscopically homogeneous 
systems, as the present case of the HS crystal. For more 
technical details, the reader is referred to references M, 
g and H in relation to TF, MWDA and GELA respec- 
tively and to reference jlj] for an extensive discussion of 
these and other functional approximations. 

The one-particle density, p(r), is usually assumed to 
be described by a sum of normalized Gaussians: 



P(r) 



(-) 

7T 



2 ^e~ Q(r ~ R) 

R 



(8) 



where R is the vector position of the crystal lattice and 
a is the Gaussian width parameter. With this density 
parametrization, the variational principle reduces to find- 
ing the value of a that minimizes the free energy func- 
tional at each mean density. Several attempts trying to 
improve the parametrization of the one-particle density 
have shown the goodness of the Gaussian one pT| . On 
the other hand, simulations have shown that deviations 
from Gaussian form are only significant at low densities 
but only at the tails of the distribution fll2| . For those 
reasons, the Gaussian parametrization seems to be an ex- 
cellent description of the one-particle density of the HS 
crystal over all range of physical densities. We shall give 
another evidence based on MC simulations |l3|-|l7|] . 

The equation of state for the HS solid is obtained as 
follows. After the minimization process we obtain the 
free energy per particle f(p) = f(p,a(p)) at each mean 
density p. Notice that a is, after the minimization, a 
function of the mean density p. Then, at fixed tempera- 
ture, the pressure is given by: 



PP R df(p) 
p dp 



(9) 



Because the standard division of the Helmholtz free en- 
ergy (0) into the ideal contribution and the excess one, 
it follows the equivalent for the free energy per particle, 
namely f = fid + fex- Therefore, the equation of state 
(0) can be formally split into two terms: 



0Pid R df ld { P ) 

= PP~. 5 

P op 



PPe, 
P 



= f3p 



d fex(p) 
dp ' 



(10) 



(11) 



where P = P L d + P ex . Following Denton et ai, we call 
Pid and P ex 'ideal-gas pressure' and 'excess pressure' re- 
spectively. Notice that Pid is not the usual ideal gas 
pressure, i.e., that which gives the ideal compressibility 
factor. In addition to the different functional approach 
used by Denton et al. for the Helmholtz free energy ex- 
cess, these authors approximate the ideal free energy per 
particle by: 



/3/«(a) = |zn(-)+3ln(A)-|, 

Z 7T Z 



(12) 



which is exact in the limit of non overlapping Gaussians. 
It is a very good approximation over all density range 
of the HS solid. However, we use the exact functional 
expression (^) to obtain the TF ideal contribution. 

An important and direct test for the equation of state 
is its capacity to predict the isothermal compressibility 
of the HS crystal. This is given by: 



pk B Txr = 



dp 



(13) 



which, extending the above formal division into an ideal 
part and an excess part, can be written as: 



pksTxid = 



pk B Txe 



dpPjd{p) 
dp 



d[3P ex {p) 
dp 



with 



11 1 

X Xid Xex 



(14) 



(15) 



(16) 



III. THE AVERAGE OF THE TWO PARTICLE 
DISTRIBUTION FUNCTION 



The averaged two-particle distribution, g(r), is defined 



as: 



WV2 



) = iJv/^ 12 / dr2p(2)(ri ' r2) ' (17) 



where V is the volume, p( 2 )(ri,r2) the two-particle den- 
sity function and df2i2 the differential solid angle aper- 
ture around T\2- In the uniform limit Eqn.([l^) reduces to 
the well known radial distribution function. MC results 
for this function were parametrized originally by Weis 
||13| with the following analytical expression: 



and 



with 



and 



g(r) =0, < r < 1 



r>l, 

i>i 



9 W (r) = ~e 



£_ [VKi(r-ri)] 2 -[VK 2 (r-n)] d 



(18) 



(19) 



(20) 



where is the number of lattice sites at distance Ri 
and A is determined by the virial theorem. The parame- 
ters ri , W\ , W2 and W are elaborated analytic functions 
of the density that give a good fitting to results of MC 
simulation (distances in all above equations are given in 
HS diameter units). These functions were refined succes- 
sively ||T^-|l6||. The parameters provided by Choi et al. 
|l7f overcome those of previous authors giving an accu- 
rate description of the MC results. The maximum root 
mean square deviation of g(r) from MC data computed 
over the distance range (r up to 3.3) is 0.17 at the highest 
density. It quickly decreases to less than 0.06 at packing 
fractions lower than r\ = 0.65. 

Here, the interesting point is the functional form used 
to fit the MC data. The function g(r) is written as a sum 
of peaks corresponding to successive shells of neighbours. 
The first peak, Eq.(|2(i|), has a characteristic form but all 
the rest, Eq. (pl|) , have exactly the same functional form. 
These later peaks only differ from each other on the dis- 
tances where they are located, Ri, and on their normal- 
izations which correspond to the number of neighbours 
at distance Ri (according to a lattice without vacancies) . 

The crucial point is to realize that if the spherical av- 
erage of p^ (r, r'), which defines g(r), is done on the 
product of two sums of Gaussians, p(r)p(r'), we obtain 
precisely all the peaks of g(r) except the first one. In 
effect, from the definition of this average: 

0o(ri2) - J dn 12 J dr 2 p( ri )p(r 2 ), (22) 

it is straightforward to obtain 

9o(r)=J29oHr), r > (23) 

i>0 



with 



9i°\r)=^-(^- )hae~^\ (24) 
47rp Zir 



and 



9o (r) = —(—)»ni 25 

4np 2tt RiT 



1 Q 1 

( — ) 2 n 4 

4-7T/9 2n 



R*r 



i > 1. 



(26) 



We have dropped the second exponential in the expres- 
sion ( p5|) as its contribution to each peak is, at most, 
20 orders of magnitude smaller than the first exponen- 
tial (for the current values of the Gaussinan parameter 
a). The first peak, Eq. (p4|) , of go(r) does not appear in 
g(r) because the exclusion of the self- interaction. The 



second one, Eq.(|26|) for i = 1, differs from the first one, 
Eg. (pp[) , of g(r). However, identifying % with W 2 , all the 
remaining peaks are exactly the same both in g(r) and 
in 9o(r). 

The identification of peaks in g(r) with those in ga(r) 
has several interesting consequences. The immediately 
obvious one is that the whole two-body correlation be- 
tween particles beyond nearest-neighbours is already in- 
cluded in the product p(r)p(r'). This point has been sug- 
gested and extensively discussed by two of us in relation 
with the perturbation weighted density approximation 
(PWDA) for simple systems with attractive interaction 
potentials Jf8|,[f9| and it is confirmed in the present dis- 
cussion. We have proposed that most of the correlation 
of the two particle distribution function, 



p( 2 \r,r') = p(r)p(r')g(v,v'), 



(27) 



is already described by the mentioned product. Thus, 
g(r, r') could be basically approximated by a step func- 
tion to exclude the self-interaction. Instead, we went fur- 
ther and mapped <?(r, r') into a homogeneous fluid at an 
very low effective density determined by the compress- 
ibility equation. This proved to be an excellent criterium 
to determine the perturbation contribution to the free 
energy @. 

On the other hand, the identification of peaks also con- 
firms the goodness of the Gaussian parametrization for 
the one-particle density. Moreover, it suggests that from 
the theoretical minimization process of the free energy 
one can obtain information on the crystal average dis- 
tribution function g(r). (A work along this line is in 
progress). Here, we explore the other way: from MC data 
(the parameter W in (|2~l|)) we obtain the corresponding 
Gaussian width parameter a of the 'MC Gaussian one- 
particle density'. Then, we are able to compute the MC 
ideal contribution to the free energy, Fid, throughout the 
exact functional form Eq. (||) . The MC free energy excess 
is now immediately obtained and so are the ideal and 
the excess contributions of the pressure and the com- 
pressibility. Notice that the excellent fitting of Choi et 
al. allows us to evaluate analytically all these thermo- 
dynamics properties over the density range up to almost 
close packing. In figure I, for a, and figure 4, for P i( i and 
P exi we show some of these results in comparison with 
those inferred by Denton et al. Q from Young and Alder 
p0| simulation data. The fair agreement is a sign of the 
consistency of the different MC data. In any case, the 
treatment of MC data we have followed is much more 
powerful as it gives a continuous expression of the ther- 
modynamic variables, as a function of density, up to al- 
most close packing. 



IV. RESULTS AND DISCUSSION 

Figure I shows the logarithm of the Gaussian parame- 
ter a inferred by us from Choi et al. MC simulations fTij ] 
and by Denton et al. M from Young and Alder MC sim- 
ulations §(|. The predictions of TF, MWDA and GEL A 
are also displayed. (MWDA data in all figures are those 
obtained by Tejero et al. which we assume free of possible 
convergence problems). TF predicts a systematic overes- 
timation of the Gaussian width parameter a but always 
in fair agreement with simulation results. GELA predic- 
tions are better at low densities as it was already known. 
At high densities both, TF and GELA, predict close val- 
ues for a. The MWDA gives the best predictions ex- 
cept at the lowest densities. Figure 2 shows the equation 
of state predicted by TF, MWDA and GELA together 
with simulation results. At low densities (see also Figure 
3), the GELA pressure is slightly above the simulation 
pressure while TF and MWDA pressures have lower val- 
ues. All of them improve their agreement with simulation 
as the pressure increases. At the highest density where 
we have performed calculations, the TF pressure differs 
from the simulation one in less than 0.3%. In Figure 4, 
we show the two contributions to the pressure, ideal and 
excess, obtained from TF, MWDA and GELA in com- 
parison with those inferred from simulation data. TF 
and GELA show an excellent agreement with simulation 
results. The MWDA results are also quite good. Notice 
that the good theoretical predictions for the a parameter 
(Figure I) will necessarily mean comparable good predic- 
tions for the ideal pressure (Figure 4) because equation 
(^) is exact. However, due to the functional approxima- 
tions proposed for the excess free energy (or for the total 
free energy) the theoretical excess pressure (or the total 
pressure) must be also compared with simulation data. 
This is even more important at high densities where a sig- 
nificant error in the excess contribution would not have 
appreciable effects on the total pressure. Figure 5 shows 
the compressibility predicted by TF together with that 
obtained from simulation. We have used the standard 
cubic spline treatment of the pressure data to obtain the 
compressibility via equation (13). Once more, the agree- 
ment is quite good. We do not have accurate enough 
data of GELA and MWDA to obtain the compressibil- 
ity properly. It exhibits unphysical oscillations which, 
most probably, are due to round off effects of pressure 
data. For completeness, we present in inset of figure 5 
the inverse of the ideal and excess compressibilities. The 
agreement with simulation is good. 

The splitting of the pressure into those parts is quite 
convenient for all the above discussion. However, we are 
inclined to think that there is no physical meaning in the 
ideal and the excess pressures, such as they are defined 
by Eq.© and (O). 



V. CONCLUSIONS 
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In this paper, we have presented two different contri- 
butions. One concerns to the method used to handle 
MC data of the g(r) and the other concerns to the accu- 
racy of the equation of state computed with functional 
approaches. 

We have proposed a new method, using MC data of 
g(r) , to determine the rms displacement of a particle from 
its equilibrium lattice site (i.e. the parameter a if p(r) 
is approximated by a sum of Gaussians). This should 
work perfectly for peaks of g(r) at distances where the 
two-particle distribution function is already given by the 
product p(y)p(y'). In practice, it is possible to deduce 
already the rms displacement (or parameter a) from the 
second peak of g(r). This procedure should be consistent 
with a direct determination of the rms displacement from 
MC configurations. Assuming the Gaussian parametriza- 
tion for p(r), an accurate fitting of g(r) allows us to deter- 
mine the parameter a as a function of the density. From 
this, the ideal pressure is straightforwardly obtained via 
the exact expression (||). From the total and the ideal 
pressures, the excess pressure follows immediately. 

As a consequence of the analytical fitting of MC data 
for g(r), it can be deduced that beyond the nearest- 
neighbours the two-particle correlation function is practi- 
cally given by the product p(r)p(r'), i.e., g(r, r') » 1. Wc 
are not aware that this quite interesting and important 
result was previously mentioned by other authors. 

The other aim of this paper has been to show the ex- 
cellent behaviour of the equation of state of the FCC 
HS solid computed with functional approaches. They 
agree with simulation results up to almost close packing. 
This remarkable behaviour is extended to the contribu- 
tions to the pressure, namely the ideal and the excess 
pressures, and to the compressibility and their ideal and 
excess parts. They all are in notable agreement with 
simulation up to almost close packing. 

Finally, we want to remark that functional approaches 
provide a reliable reference HS system for perturba- 
tion theories, especially at high densities. This makes 
them particularly suitable for describing the solid-solid 
isostructural transition of simple systems with extremely 
short ranged attractive potentials [0]. We have applied 
the PWDA Jl8| , p^| mentioned above to an attractive 
square well |2l| ] and to HS plus Yukawa attractive tail 
p2fl . Precisely, the PWDA is based, among other things, 
on two of the points explored here: the goodness of the 
TF, even at high densities, for describing HS and the ac- 
curacy of the product p(r)p(r') for describing the correla- 
tion between particles beyond nearest neighbours. Note 
that after Tejero et al. discussion on the numerical prob- 
lems at high densities of the MWDA, results based on 
this functional approach must be seem with some cau- 
tion g|. 
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FIG. 1. Logarithm of the Gaussian width parameter a vs. 
packing fraction r\. Solid curve inferred, in this paper, from 
simulation data of Choi et al. Solid squares inferred by Den- 
ton et al. from simulation data of Young and Alder. The 
dashed curve is the TF prediction. Solid triangles and open 
triangles are the results from GELA and MWDA respectively 
obtained by Tejero et al.. 



FIG. 2. Logarithm of the equation of state Pj ' pksT vs. 
packing fraction rj. Curves and symbols as in figure 1. 

FIG. 3. Equation of state P/pk B T vs. packing fraction rj 
at low densities. Curves and symbols as in figure 1. 

FIG. 4. Equation of state P/pksT vs. packing fraction rj. 
Curves and symbols as in figure 1. 

FIG. 5. pksTxr vs. packing fraction rj. Inset: 1/ ' pkaTxr 
vs. packing fraction rj. Curves as in figure 1. 



